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ABSTRACT 


The  sensitivity  of  numerical  solutions  of  systems  of  nonlinear  flow  equations 
(Navier-Stokes  equations)  to  the  grid  used  is  investigated  through  the  use  of  the  South 
China  Sea  (SCS)  numerical  ocean  model.  Traditionally,  rectangular  coordinate  grids  are 
used  in  environmental  modeling.  The  advantage  of  rectangular  coordinate  grids  is  their 
simplicity  in  the  generation  process.  However,  rectangular  coordinate  grids  are  not  well 
suited  for  regions  with  complex  terrain  (coastlines  and  topography)  and  occasionally  lead 
to  poor  accuracy  in  numerical  solutions.  The  grid  generation  techniques  are  being 
introduced  to  coastal  ocean  modeling  to  study  the  sensitivity  of  numerical  solutions  to  the 
grid  used  and  to  investigate  the  enhancement  of  the  modeling  process.  Grid  generation 
techniques  are  broadly  used  in  the  aeronautical  engineering  community  for  solving  CFD 
problems. 

One  orthogonal  (121x191)  and  two  curvilinear  nearly-orthogonal  grids  (121x191 
and  151x241)  are  designed  to  couple  with  the  SCS  numerical  ocean  model.  The  grids  are 
designed  using  the  EAGLEView  grid  generation  code  developed  by  the  National  Science 
Foundation  (NSF)  Engineering  Research  Center  (ERC)  of  Mississippi  State  University. 
EAGLEView  implements  a  grid  generation  technique  using  mainly  elliptic  and  algebraic 
generation  systems.  The  designed  grids  are  processed  with  the  SCS  numerical  ocean  model 
for  200  days  to  study  the  sensitivity  of  numerical  solutions  to  the  grid  used.  The  solutions 
of  the  temperature  and  salinity  fields  are  presented  and  analyzed.  The  advantages  of 
curvilinear  nearly-orthogonal  grids  are  also  discussed. 
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I.  INTRODUCTION 


A.  GENERAL  REMARKS 

The  advent  of  supercomputers  with  high  computing  speed  and  large  core 
memories  has  allowed  the  majority  of  unsolved  problems  in  Fluid  Dynamics  to  be  treated 
numerically.  Complex  problems  can  now  be  solved  in  a  few  seconds  of  computer  time  with 
very  little  cost  which  would  have  taken  years  of  work  with  computational  methods  and 
computers  available  twenty  years  ago.  As  a  consequence,  a  new  method  termed 
computational  fluid  dynamics  (CFD)  has  become  available  for  attacking  complicated 
problems  in  systems  of  non-linear  partial  differential  equations.  This  method  has 
permitted  advanced  simulations  of  flow  phenomena  for  a  variety  of  applications.  The 
areas  range  from  aircraft  and  missile  design  to  large-scale  simulations  of  the  atmosphere 
and  the  ocean.  The  numerical  solutions  of  the  fluid  flow  equations  have  progressed  to  the 
point  of  finite-difference  (or  finite-volume)  approximation  about  arbitrary  and  complex 
geometries  with  reliable  and  predictable  accuracy.  Similar  advances  are  also  emerging  in 
hydrodynamics,  electromagnetism,  magnetohydrodynamics,  heat  and  mass  transfer  and,  to 
some  degree,  in  all  field  problems. 

The  traditional  experimental  and  theoretical  approaches  remain  important  in 
engineering  design  for  systems  involving  fluid  and  heat  flow  equations.  However,  the 
trend  is  clearly  toward  greater  reliance  on  computer  based  predictions  in  the  design  process. 
This  trend  can  be  largely  explained  by  the  effect  of  a  rapidly  developing  digital  computer 
industry  and  a  relatively  low  cost  of  utilizing  these  machines.  The  cost  of  performing  a 
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given  calculation  has  been  significantly  reduced  in  the  last  decade  and  is  expected  to  be 
steadily  declining  for  some  time  into  the  future. 

Contrarily,  the  experimentation  costs  have  been  significantly  increased  in  recent 
years.  The  construction  of  complex  models  and  wind  tunnel  facilities  to  match  realistic 
flow  conditions  has  elevated  the  cost  every  day.  The  increasing  energy  cost  to  operate  the 
wind  tunnels  is  also  contributing  to  the  mounting  expense  of  the  experimental  approach. 
However,  the  experimental  approach  still  remains  the  best  method  of  representing  the  most 
realistic  solutions  for  many  complex  flow  problems. 

On  the  other  hand,  the  theoretical  approach  is  the  most  inexpensive  method  to  solve 
flow  problems.  The  difficulty  with  the  theoretical  approach  is  that  simplifying  assumptions 
are  made  and  the  problems  are  restricted  to  simple  geometries  and  physics.  One  significant 
advantage  of  the  theoretical  approach  is  that  reasonable  answers  can  be  obtained  in  a 
minimum  amount  of  time  and  this  is  quite  useful  in  preliminary  design  work.  In  closing, 
the  computational  approach  is  judged  to  be  the  most  efficient  method  for  treating  systems 
of  non-linear  flow  equations  (Navier-Stokes  equations)  at  the  present  time. 

B.  GRID  GENERATION 

A  vital  element  of  computing  numerical  solutions  for  systems  of  non-linear  flow 
equations  is  the  discretization  of  the  physical  field  into  a  collection  of  points  or  elemental 
volumes  (cells)  on  which  to  represent  the  flow  equations  in  finite  form.  The  application  of 
finite  mathematical  techniques  to  acquire  this  essential  measure  is  called  grid  generation. 
A  grid  which  is  not  well  suited  to  the  problem  can  lead  to  unsatisfactory  results.  In  some 
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applications,  improper  choice  of  grid  point  locations  can  also  lead  to  model  instability  or 
lack  of  convergence  in  numerical  solutions  [1],  A  well  constructed  grid  can  greatly 
simplify  solutions  of  systems  of  non-linear  fluid  flow  equations  and  subsequently  reduce 
computing  time  and  the  cost  of  performing  the  calculation. 

Grid  generation  is  a  construction  process  which  must  accommodate  the 
characteristic  of  the  governing  equations  and  the  geometry  of  the  physical  domain.  The 
governing  equations  are  defined  in  a  physical  coordinate  system  and  are  then  transformed 
to  an  idealized  rectangular  computational  domain.  This  situation  is  illustrated  in  Figure  1. 
In  simple  terms,  grid  generation  is  a  procedure  for  the  orderly  distribution  of  observers  over 
a  physical  field  in  such  a  way  that  efficient  communication  between  observers  is  possible 
and  all  physical  phenomena  on  the  entire  continuous  field  may  be  represented  with 
sufficient  accuracy  by  this  finite  collection  of  observations  [2].  Hence,  the  representation 
of  the  physical  phenomena  on  the  continuous  field  is  strongly  dependent  on  the  distribution 
of  these  observers.  This  implies  that  numerical  solutions  of  nonlinear  flow  equations  is 
heavily  affected  by  the  distribution  of  the  grid  points. 

Grid  generation  about  complex  configurations  is  considered  as  a  major  pacing 
factor  for  the  ability  to  obtain  numerical  solutions  of  non-linear  partial  differential 
equations.  There  are  generally  three  available  methods  to  generate  numerical  grids:  (1) 
complex  variable  method,  (2)  algebraic  method,  and  (3)  differential  equation  method  based 
on  the  solution  of  a  set  of  partial  differential  equations  [1],  The  first  method  has  the 
advantage  of  analyzing  or  partially  analyzing  the  transformations  that  are  being  used  as 
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Physical  Domain  Computational  Domain 


Figure  1:  Physical  domain  and  computational  domain  [3] 


opposed  to  the  other  two  methods  that  are  entirely  numerical.  Unfortunately,  this  method 
is  restricted  solely  to  two  dimensions  and  therefore,  it  is  uncommonly  used.  The  second 
method  is  based  on  the  interpolation  or  approximation  among  boundaries  and  intermediate 
surfaces  in  the  field.  The  last  method,  in  which  the  coordinates  are  the  solution  of  partial 
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differential  equations,  may  be  elliptic,  parabolic,  or  hyperbolic.  Several  major  grid 
generation  codes,  such  as  EAGLEView,  GENIE  ++,  INGRID,  GRAPE,  and  GRIDGEN, 
are  presently  available  for  use  to  generate  grids  on  complex  configurations.  The  techniques 
of  numerical  grid  generation  and  its  applications  are  thoroughly  covered  in  reference  [4], 

C.  THE  STUDY  OF  NUMERICAL  SOLUTIONS 

The  intent  of  this  thesis  is  to  study  the  sensitivity  of  numerical  solutions  of  systems 
of  nonlinear  flow  equations  due  to  grid  utilization.  The  study  involved  the  design  of  an 
orthogonal  and  two  curvilinear  nearly-orthogonal  grids  which  were  used  in  conjunction 
with  the  South  China  Sea  (SCS)  numerical  ocean  model.  The  grids  were  designed  using 
the  EAGLEView  code  which  was  developed  by  the  National  Science  Foundation  (NSF) 
Engineering  Research  Center  (ERC)  for  Computational  Field  Simulation  of  Mississippi 
State  University.  The  SCS  numerical  ocean  model  was  developed  at  the  Naval 
Postgraduate  School  (NPS)  based  on  the  Princeton  Ocean  Model  (POM)  [5],  The 
sensitivity  of  numerical  solutions  of  the  SCS  model  due  to  the  design  of  the  three  grids  is 
the  scope  of  this  study.  The  advantages  and  disadvantages  of  orthogonal  and  curvilinear 
grids  will  also  be  discussed. 

Traditionally,  rectangular  coordinate  grids  have  been  commonly  used  in 
environmental  modeling.  The  advantage  of  rectangular  coordinate  grids  is  their  simplicity 
in  the  generation  process.  However,  open  boundary  conditions  are  poorly  represented 
numerically  with  rectangular  coordinate  grids.  This  is  especially  true  for  regions  with 
complex  terrain  even  when  treated  with  very  high  resolution.  Rectangular  coordinate  grids 
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also  lack  the  ability  to  follow  the  coastlines  and  eliminate  apathetic  regions  of  the  domain. 
This  has  resulted  in  an  unnecessary  expense  of  computing  resources  in  regions  of  apathy 
and  inaccurate  numerical  solutions  in  regions  of  interest.  On  the  other  hand,  the  curvilinear 
nearly-orthogonal  grids  have  the  ability  to  enhance  the  numerical  solutions  by  better 
treatment  of  coastlines  and  boundary  conditions  and  still  remain  within  realistic 
computational  resources.  The  curvilinear  grids  also  have  the  ability  to  increase  resolution 
in  the  subregion  of  the  model  domain  without  increasing  computational  expenses.  The  grid 
generation  technique  has  been  widely  used  in  the  aerospace  community  and  the 
introduction  of  its  application  to  the  oceanography  community  has  created  a  new  direction 
in  coastal  ocean  modeling  [6,7]. 
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II.  GRID  GENERATION  SYSTEM-THE  EAGLEVIEW  CODE 
A.  GENERAL  OVERVIEW 

The  EAGLE  View  grid  code  is  a  general  three-dimensional  grid  generation  code 
based  on  a  composite  block  structure.  The  code  was  developed  by  the  NSF  Engineering 
Research  Center  for  Computational  Field  Simulation  of  Mississippi  State  University.  The 
code  can  be  operated  either  as  an  algebraic  generation  system  or  as  an  elliptic  generation 
system.  The  code  can  also  be  operated  in  two  or  three  dimensions  on  a  plane  or  on  a  curved 
surface.  In  the  case  of  a  curved  surface,  the  surface  is  splined  and  the  generation  is  done  in 
terms  of  surface  parametric  coordinates.  The  input  of  the  code  is  structured  to  be  user- 
friendly  oriented  and  arbitrary  block  configurations  can  be  treated.  The  EAGLEView  code 
was  developed  with  the  intention  of  reducing  the  time  expended  to  generate  the  grids  for 
engineering  designs  and  analysis  in  computational  field  simulation. 

The  EAGLEView  grid  code  permits  an  arbitrary  three-dimensional  region  to  be 
filled  with  a  number  of  blocks.  The  number  of  blocks  can  be  varied  depending  on  the 
discretion  of  the  user.  The  blocks  are  linked  to  each  other  with  continuity  across  the 
interfaces.  The  code  uses  an  elliptic  generation  system  with  automatic  evaluation  of  the 
control  functions.  The  evaluation  of  the  control  functions  is  completed  either  directly  from 
the  initial  algebraic  grid  and  then  smoothed,  or  by  interpolation  from  the  boundary-point 
distributions.  In  the  former  case,  the  relative  spacing  of  the  grids  is  done  by  the  algebraic 
system  and  the  smoothing  of  the  grid  is  done  by  the  elliptic  system.  In  the  latter  case,  the 
arc  length  and  the  curvature  contributions  to  the  control  functions  are  evaluated  and 
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interpolated  separately  into  the  field  from  the  appropriate  boundaries.  The  control  function 
at  each  point  in  the  field  is  then  formed  by  combining  the  interpolated  elements. 

The  control  functions  can  also  be  automatically  determined  to  provide 
orthogonality  at  boundaries  with  specified  normal  spacing.  In  this  case,  the  iterative 
adjustments  in  the  control  functions  are  made  by  increments  radiated  from  boundary  points 
where  orthogonality  has  not  yet  been  attained.  This  allows  the  basic  control  function 
structure  evaluated  from  the  algebraic  grid,  or  from  the  boundary-point  distributions,  to  be 
retained  and  thus  relieves  the  iterative  process  from  the  need  to  establish  this  basic  form  of 
the  control  functions. 

The  generation  of  the  initial  solution  to  start  the  iterative  process  of  the  elliptic 
generation  system  is  done  by  the  three-dimensional  algebraic  system  which  is  based  on 
transfinite  interpolation  (using  either  Lagrange  or  Hermite  interpolation)  [4],  This  feature 
permits  the  code  to  be  run  as  an  algebraic  generation  system  if  desired.  The  interpolation 
is  defaulted  to  complete  transfinite  interpolation  from  all  boundaries;  however,  it  can  be 
restricted  by  input  to  any  combination  of  directions  or  lesser  degree  of  interpolation.  The 
mathematical  techniques  of  grid  generation  that  are  applicable  to  the  EAGLEView  grid 
code  are  discussed  in  detail  in  references  [4,8]. 

B.  ELLIPTIC  GENERATION  SYSTEM 

The  properties  of  the  elliptic  generation  system  is  based  on  the  derivation  of  the 
Laplace  or  Poisson  equations.  The  Laplace  system  is  the  most  simple  elliptic  partial 
differential  system  and  one  that  produces  the  smoothest  possible  grid;  however,  the 
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coordinate  lines  in  the  interior  of  the  field  tend  to  be  equally  spaced  regardless  of  how  the 
boundary  points  are  spaced.  This  is  caused  by  the  strong  smoothing  effect  of  the  Laplacian 
system.  On  the  other  hand,  the  Poisson  system  can  provide  control  of  the  coordinate  line 
spacing  in  the  field  to  match  with  that  of  the  boundaries. 

The  elliptic  grid  generation  technique  used  in  the  EAGLEView  grid  code  is  derived 
from  the  Poisson  equations 

=  0  -  1,2,3)  E,1 

where  the  functions,  P‘  are  the  “control  functions”  which  can  be  fashioned  to  control  the 
spacing  and  orientation  of  coordinate  lines. 


Reference  [4]  reveals  that  if  a  curvilinear  coordinate  system,  ^  =  2’  3  )  , 

which  satisfies  the  Laplace  system,  V2^  =  0  (z  =  1, 2,  3  )  ,  is  transformed  to  another 

coordinate  system,  (z  =  1 ,  2, 3  )  ,  then  the  new  curvilinear  coordinates,  ,  will  satisfy 
the  inhomogeneous  elliptic  system  as  defined  by  (Eq  1)  with  the  control  functions 

Pl  =  Yj  YdkP)k  0  =  1,2,3)  Eq  2 

j  =  Ik  =  1 

where  defined  the  transformation  from  4*  to  ^ , 


and 


y  y  a4md4"aV 

m=xn=xd^d^md^ 


Eq  3 


/  =  ^-nk  =  hgjmghn-gjnghm) 
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with  (i,j,h),  (k,m,n)  cyclic,  is  the  contravariant  metric  tensor. 
g  =  det\gtj\  =  •  (}_,  x  £  3  )  is  the  square  of  the  Jacobian  of  the  transformation,  and 

gy  =  £  /  are  the  elements  of  the  covariant  metric  tensor.  In  these  relations, 

^  =  xi+yj  +  zk  is  the  position  vector  of  the  grid  point  in  the  Cartesian  coordinates,  and 

4*  (i  =  1,  2,  3  )  are  the  three  new  curvilinear  coordinates. 

These  results  reveal  that  a  coordinate  system  obtained  by  subsequent 
transformation  (often  called  “stretching”  transformation)  from  one  generated  as  the 
solution  of  the  Laplace  system  can  be  solved  directly  as  the  solution  of  the  Poisson  system, 

(Eq  1),  with  appropriate  control  functions,  as  defined  by  (Eq  3).  Therefore,  the 

Poisson  system  can  be  adopted  as  the  appropriate  generation  system  with  the  control 
functions  specified  directly  rather  than  through  a  subsequent  transformation.  Thus,  the 
appropriate  generation  system  can  be  defined  by  the  combination  of  (Eq  1)  and  (Eq  2) 

(i  -  1.2,3)  Eq  5 

j  =  1  *  =  1 

The  basis  of  (Eq  5)  is  that  it  produces  a  coordinate  system  that  corresponds  to  the 
subsequent  application  of  a  stretching  transformation  from  a  coordinate  system  generated 

for  maximum  smoothness.  In  (Eq  3),  the  three  control  functions,  Pln  (i  =  1, 2,  3  )  , 
correspond  to  one-dimensional  stretching  in  each  coordinate  direction  and  are  the  most 
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important  of  the  control  functions.  The  other  control  functions  are  taken  to  be  zero,  i.e. 
pjk  =  WV  Thus,  the  generation  system  becomes 

=  g"p.  (i  =  1,2,3  )  Eq  6 

Further  transformation  of  (Eq  6)  yields 

2  Z  +  Z  gkkp/}p =  0  Eq  7 

*  =  iy  =  i  /t=i 

which  is  the  elliptic  generation  system  equation  that  is  used  in  the  EAGLEView  grid  code. 
The  control  functions  P k  are  used  to  control  the  spacing  and  orientation  of  gridlines  in  the 

interior  of  the  field.  Procedures  for  the  determination  of  these  control  functions  such  that 
(Eq  7)  reflects  the  spacing  of  an  algebraic  grid  or  the  boundary  point  spacing  in  the  field 
are  specifically  discussed  in  reference  [8]. 
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III.  SCS  NUMERICAL  OCEAN  MODEL 


A.  BASIC  MODEL  EQUATIONS 


The  South  China  Sea  (SCS)  numerical  ocean  model  was  developed  at  the  Naval 
Postgraduate  School  (NPS)  based  on  the  Princeton  Ocean  Model  (POM)  [5].  The  POM 
was  originally  created  at  Princeton  University  in  1977  for  application  to  oceanographic 
problems.  Many  subsequent  contributions  to  the  model  were  made  throughout  the  years  by 
various  scientists  and  researchers  from  the  Atmospheric  and  Oceanic  Sciences  community 
[9]. 

The  model  is  a  three-dimensional,  primitive  equation,  ocean  circulation  numerical 
model  with  realistic  coastlines,  topography  and  complete  thermodynamics.  The  model  also 
includes  an  upper  free  surface  and  16  vertical  levels  of  sigma  coordinate.  A  detailed 
description  of  the  model  is  given  in  reference  [10].  The  sigma  coordinate  system  is 
illustrated  in  Figure  2  and  is  defined  by  the  equation 


(x,  v,  t 


H(x,y) +r\(x,y,t) 

where  H(x,y  )  is  the  bottom  topography  and  q  (x,.y,  t  )  is  the  free  surface  elevation. 


Thus,  cr  ranges  from  a  =  0  at  z  =  q  to  a  =  -1  at  z  =  -H .  The  coordinates 


(x,  y,z)  represent  the  conventional  cartesian  coordinates  with  x  in  the  easterly  direction, 
y  in  the  northerly  direction  and  z  is  normal  to  the  surface  pointing  upward. 

References  [11,12,13]  indicate  that  the  velocity,  surface  elevation,  salinity  and 
temperature  fields  of  the  ocean  are  described  by  the  following  equations 
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po  :  reference  density 
g  :  gravitational  acceleration 

P  =  PeO^S)  :  in  situ  density  in  general  form  computed  by  using  the  equation  of 
state 

2Q  x  V  :  Coriolis  force  with  Q  as  the  earth  rotation  vector 

Km,  Kh  :  vertical  turbulence  exchange  coefficients  for  momentum,  heat  &  salt 

0;  :  mean  potential  temperature 

F,  Fq  :  horizontal  mixing  terms 

The  horizontal  mixing  terms  F  (Fx,  Fy  )  in  (Eq  9)  and  F0  in  (Eq  12)  can  be  written 
as  follow 
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turbulence  closure  scheme.  The  turbulent  closure  scheme  is  characterized  by  the  turbulent 

2 

kinetic  energy  (TKE),  q  /  2,  and  the  turbulent  mixing  length,  l.  Thus,  both  of  the 
turbulent  equations  can  be  written  in  the  same  form  by  the  following  equation 


KmQ  =  Km,Kkg  =  IE,Kk;Ke  -  W-Fg  =  F,  Eq  18 

The  term  W  in  (Eq  18)  is  the  “wall  proximity”  function.  In  (Eq  16),  the  first  term  on  the 
right  hand  side  represents  diffusion,  the  second  two  terms  are  shear  production,  followed 
by  buoyancy,  dissipation,  and  the  last  term  models  horizontal  mixing  [11,13].  The  vertical 
turbulent  exchange  coefficients  (TEC)  Km,  Kh,  Kq  in  the  above  equations  are  defined  as 
follows 

Km  =  l(lSm>Kh  =  1(lSh’Kq  =  Wq  Eq  19 

where  Sm,  Sh’Sg  are  the  stability  functions  and  are  given  in  reference  [14]. 
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The  model  assumed  that  the  ocean  is  hydrostatic  and  incompressible  (Boussinesq 
approximation).  Wind  stress,  heat,  and  salinity  fluxes  are  prescribed  at  the  free  surface 
where  z  =  r\  (x,y  )  .  Zero  heat  and  salinity  fluxes  are  designated  at  the  bottom  where 
z  =  -H  (x,  y  )  .  The  condition  of  no  diffusive  fluxes  of  any  property  across  the  interface 
is  used  at  land  boundaries.  The  condition  of  mode  splitting  of  the  external  (barotropic)  and 
internal  (baroclinic)  modes  is  imposed  in  the  model  to  enhance  computer  resources  [11,13]. 

B.  MODEL  INITIALIZATION 

The  model  is  initialized  with  10  seconds  time  step  for  the  external  (barotropic) 
mode  and  400  seconds  time  step  for  the  internal  (baroclinic)  mode  to  satisfy  the  Courant- 
Friedrichs-Levy  (CFL)  computational  stability  criterion.  The  bathymetry  is  obtained  from 
the  global  ET0P05  bathymetry  dataset  with  5’  x  5’  resolution  and  is  interpolated  to  the 
grid  of  the  model.  The  bathymetry  is  filtered  to  remove  high  frequency  noise.  The  SCS 
model  is  simulated  with  close  boundaries  in  order  to  simplify  the  modeling  process  and  to 
strictly  focus  on  the  sensitivity  study  of  numerical  solutions  due  to  grid  usage.  The  model 
is  also  initialized  with  Levitus  salinity  and  temperature  climatology  [15].  The  annual 
salinity  and  temperature  fields  from  the  Levitus  dataset  are  interpolated  into  the  model 
grids  and  spun  up  30  days  diagnostically  (salinity  and  temperature  fields  held  constant) 
prior  to  making  the  prognostic  run. 
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IV.  NUMERICAL  SIMULATIONS 


A.  GRID  DESIGN 

One  orthogonal  grid  (gridl)  and  two  curvilinear  nearly-orthogonal  grids  (grid2  and 
grid3)  are  designed  to  cover  the  entire  domain  of  SCS.  Gridl  (Figure  3)  and  grid2  (Figure 
4)  are  created  to  have  the  same  number  of  grid  points  (121  x  191  =  23,1 1 1  points).  Grid3 
(Figure  5)  is  designed  to  have  a  much  higher  resolution  of  grid  points  (151  x  241  =  36,391 
points).  The  traditional  rectangular  grid,  gridl,  is  completely  overlaid  on  the  SCS  domain 
from  99°  -  121  °East  longitude  and  from  0°  -  25° North  latitude.  This  produces  a  range 
resolution  of  approximately  20  km  x  15  km  for  gridl.  Grid2  and  grid3  are  designed  to 
eliminate  apathetic  regions  of  the  SCS  domain  as  much  as  possible  without  causing  the 
effect  of  skewness  to  the  grids  (grids  are  ineffective  to  model  the  flow  equations  if  the  level 
of  skewness  is  too  high).  The  average  range  resolution  for  grid2  is  approximately  12  km  x 
15  km  and  for  grid3  is  approximately  10  km  x  12  km.  The  subregions  near  the  northern 
boundary,  25° North  latitude,  and  the  southern  boundary,  0°  latitude  (the  equator),  of 
grid2  and  grid3  have  a  higher  horizontal  resolution  from  east  to  west  than  the  central  part 
of  the  domain.  The  average  range  resolution  for  these  subregions  is  approximately  7  km  x 
1 5  km  for  grid2  and  5  km  x  1 2  km  for  grid3 .  This  is  one  advantage  of  the  curvilinear  nearly- 
orthogonal  grids,  which  allows  the  resolution  of  the  subregions  of  the  domain  to  be 
increased  without  increasing  the  matrix  size  of  the  whole  domain. 

The  above  designed  grids  are  processed  with  the  SCS  numerical  ocean  model  for 
200  days  to  analyze  the  sensitivity  of  numerical  solutions  due  to  grid  used.  The  temperature 
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and  salinity  fields  of  the  model  are  initialized  by  reference  [15].  The  model  is  spun  up  for 
30  days  to  obtain  the  initial  solutions  of  temperature  and  salinity  fields  (0  days  simulations). 
Subsequent  solutions  are  obtained  for  30  days,  50  days,  100  days,  150  days  and  200  days 
respectively.  The  simulations  are  computed  on  the  NPS  CRAY-YMP  supercomputer  using 
multitasking  modes  with  8  CPU  operated  correspondingly.  The  numerical  solutions  of  the 
higher  resolution  grid,  grid3,  are  used  as  a  reference  for  the  true  solutions.  The  numerical 
solutions  of  the  surface  temperature  (sigma  level  1)  and  salinity  fields  of  the  three  grids 
(gridl,  grid2  and  grid3)  are  compared  and  analyzed.  The  solutions  of  sigma  level  10  and 
level  15  are  also  obtained  and  compared  for  all  three  grids.  The  solutions  of  grid2  are 
expected  to  be  closely  related  to  the  true  solutions  (i.e.  solutions  of  grid3). 

B.  ANALYSIS  OF  NUMERICAL  SOLUTIONS 

The  numerical  solutions  of  the  salinity  and  temperature  fields  are  visualized  by 
using  the  Flow  Analysis  Software  Toolkit  (FAST).  FAST  is  a  graphic  software  package 
developed  at  NASA  Ames  Research  Center,  Moffett  Field,  CA  and  is  specifically  designed 
for  scientific  visualization  in  computational  field  simulation.  In  the  near  future,  a  modified 
version  of  EAGLEView  will  be  included  as  an  integrated  module  in  FAST  to  provide  the 
researcher  with  a  fully  integrated  system  from  the  design  process  to  system  simulation  and 
visualization  [17,18,19,20],  The  visualizations  of  the  salinity  and  temperature  fields  for 
sigma  level  1,10  and  15  of  the  three  grids  are  shown  in  Appendix  A  and  B.  Appendix  A 
contains  the  simulations  of  the  salinity  field  and  Appendix  B  contains  the  simulations  of 
the  temperature  field  for  0  days,  100  days  and  200  days  respectively. 
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1.  Salinity  Field 

The  numerical  solutions  of  the  three  grids  (gridl,  grid2  and  grid3)  for  sigma  level 
1  (surface  level),  level  10  and  level  15  are  plotted  to  allow  the  visualization  of  the  salinity 
field.  The  initial  solutions  (0  days  simulations)  of  the  surface  level  for  all  three  grids  are 
closely  related  (Figures  Al,  A2,  A3).  However,  there  are  distinct  differences  in  the  flow 
pattern  between  the  200  day  simulations  of  the  three  grids  (Figures  A4,  A5,  A6).  The 
differences  are  more  noticeable  in  the  subregions  near  the  northern  boundary  of  the  SCS 
domain  where  the  salt  density  is  high  (~  34.1  part  per  thousand  -  ppt)  compared  to  other 
parts  of  the  domain.  The  solution  of  grid2  is  likely  to  resemble  the  solution  of  grid3 
(Figures  A5,  A6).  Contrarily,  the  solution  of  the  traditional  rectangular  gridl(Figure  A4), 
is  distinguished  from  the  solutions  of  grid2  and  grid3.  The  flow  pattern  of  high  salt  density 
(~  34.1  ppt)  of  gridl  falls  short  of  the  Hainan  island  while  the  same  flow  pattern  of  grid2 
and  grid3  extends  beyond  the  island  into  the  Gulf  of  Tonkin.  The  salinity  flow  patterns  of 
the  SCS  central  and  southern  region  of  the  domain  are  also  different.  The  magnifications 
of  the  northern  region  of  the  SCS  domain  for  0  days,  100  days  and  200  days  simulations 
are  presented  in  Figures  A7  to  A15.  Notice  that  the  flow  pattern  of  grid2  for  the  100  day 
simulation  also  more  closely  resembles  the  true  solution  (solution  of  grid3)  compared  to  the 
solution  of  gridl  (Figures  A 10,  Al  1,  A 12). 

The  differences  in  these  flow  patterns  between  the  solutions  of  the  three  grids  are 
caused  by  the  differences  in  the  spatial  distribution  of  the  grid  points  (i.e.  Ax,  Ay).  The 
curvilinear  grids  (grid2  and  grid3)  have  a  higher  range  resolutions  (Ax,  Ay)  in  regions  of 
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interest  than  the  traditional  rectangular  gridl  and  this  subsequently  causes  the  sensitivities 
in  numerical  solutions  between  the  three  grids  as  discussed  above. 

The  200  day  simulations  of  level  10  (Figures  A16  to  A18)  and  level  15  (Figures 
A19  to  A21)  are  also  magnified  in  the  northwestern  region  of  the  domain  to  analyze  the 
sensitivity  of  the  salinity  field  below  the  sea  surface.  The  solutions  of  gridl  are  the  most 
different  from  the  solutions  of  grid2  and  grid3  for  both  level  10  and  level  15  as  expected 
since  gridl  has  a  much  smaller  range  resolution  than  grid2  and  grid3.  The  density  of 
salinity  along  the  coast  of  southeastern  China  is  smaller  for  gridl(Figures  A16,  A19) 
compared  to  grid2  (Figures  A 17,  A20)  and  grid3  (Figures  A 18,  A21).  The  flow  patterns 
inside  the  Gulf  of  Tonkin  similar  between  grid2  and  grid3.  The  density  of  the  salinity  field 
for  this  same  flow  pattern  of  gridl  is  also  smaller  than  for  grid2  and  grid3.  These 
sensitivities  in  numerical  solutions  can  also  be  explained  by  the  differences  in  grid  point 
distributions  between  the  three  grids  as  in  the  case  of  the  surface  level. 

2.  Temperature  Field 

The  solutions  of  the  temperature  field  are  similarly  plotted  for  all  three  sigma  levels 
(level  1, 10  and  15).  The  initial  surface  temperatures  (0  days  simulations)  for  all  three  grids 
are  very  closely  related  as  in  the  case  of  the  salinity  field  (Figures  Bl,  B2,  B3).  The 
simulations  of  the  200  day  surface  temperature  field  clearly  show  cold  water  mass  (~ 
22 0  C)  flowing  from  the  north  into  SCS  (Figures  B4,  B5,  B6).  Similar  to  the  salinity  field, 
the  flow  pattern  of  these  cold  water  masses  extends  deep  into  the  domain  of  the  SCS  from 
east  to  west  for  grid2  and  grid3  and  falls  short  for  gridl .  The  northern  region  of  the  domain 
is  magnified  for  the  solutions  of  0  days,  100  days  and  200  days  respectively  (Figures  B7  to 
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B 1 5).  The  flow  pattern  of  the  cold  water  mass  (~  22  0  C)  for  the  1 00  day  simulation  of  grid2 
is  more  closely  related  to  the  simulation  of  grid3.  The  flow  of  the  cold  water  mass  of  gridl 
also  falls  short  for  the  100  day  simulation  similar  to  the  200  day  simulation  as  discussed 
above. 

The  200  day  simulations  of  the  northeastern  region  of  the  SCS  domain  is  magnified 
for  level  10  (Figures  B16  to  B18)  and  level  15  (Figures  B19  to  B21)  to  investigate  the 
temperature  field  deep  below  the  sea  surface  level.  The  features  of  the  temperature  field  at 
these  two  levels  are  difficult  to  distinguish.  However,  cold  water  (~  16°  C  for  level  10  and 
~  12°  C  for  level  15)  appears  at  the  northern  boundary  between  China  and  Taiwan  for 
solutions  of  gridl ,  where  the  temperature  at  this  region  for  grid2  and  grid3  is  ~  22°  C  for 
level  10  and  20°  C  for  level  15.  Similar  to  the  salinity  field,  the  solutions  of  grid2  are  more 
closely  related  to  the  reference  solutions  of  grid3. 
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y.  DISCUSSION  AND  CONCLUSION 


The  numerical  solutions  of  the  system  of  non-linear  flow  equations  (Navier-Stokes 
equations)  are  very  sensitive  to  grid  designs.  This  is  notably  true  in  coastal  ocean  modeling 
where  large  domains  of  the  ocean  are  simulated.  Poorly  designed  grids  can  lead  to 
inaccurate  representations  of  the  physical  field  and  an  unnecessary  expense  of  computing 
resources.  Gridl  and  grid2  have  the  same  number  of  grid  points  (121  x  191  =  23,111 
points);  however,  the  rectangular  gridl  has  only  13,005  effective  grid  points  (i.e.  the 
number  of  grid  points  over  the  water  region  of  the  domain),  56%  of  23,1 1  lpoints,  that  are 

useful  in  the  calculation  of  numerical  solutions.  Grid2  has  16,252  effective  grid  points, 

\ 

70%  of  23,1 1 1  points,  that  are  used  in  the  simulation  of  the  fluid  flow.  Notice  that  a  certain 
part  of  the  apathetic  water  regions  has  been  eliminated  from  gridl  when  grid2  is  designed 
(i.e.  the  water  regions  of  the  Strait  of  Malacca  and  the  Zulu  Sea).  This  implies  that  the 
number  of  effective  grid  points  for  the  rectangular  gridl  (Figure  3)  is  much  less  than  1 3,005 
points  («  56%)  for  the  region  of  interest  (i.e.  the  SCS  domain  as  shown  in  the  curvilinear 
grid2,  Figure  4).  The  curvilinear  nearly-orthogonal  grid  has  more  useful  grid  points  for  the 
simulation  than  the  traditional  rectangular  grid  of  the  same  matrix  size.  The  curvilinear 
grid  also  has  the  capability  to  enhance  the  numerical  solutions  by  increasing  the  horizontal 
resolution  of  the  interested  subregions  of  the  domain  without  increasing  the  matrix  size  of 
the  grid  points  (the  subregions  near  the  northern  and  southern  boundary  of  the  SCS  domain 
have  a  much  higher  horizontal  range  resolution  compared  to  the  same  regions  of  gridl). 
The  CPU  time  for  grid2  is  slightly  higher  than  gridl;  however,  numerical  solutions  are 
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greatly  enhanced  as  compared  to  the  solutions  of  a  much  higher  resolution  grid3.  This  also 
directly  results  in  the  achievement  of  reduced  computational  expenses. 

The  study  produces  some  important  elements  in  numerical  simulations,  especially  in 
coastal  ocean  modeling.  The  interpretation  of  the  physical  fields  is  strongly  affected  by 
grid  point  distributions  (i.e.  grid  resolutions).  The  study  also  reveals  that  the  grid 
generation  technique  has  some  advantages  over  the  traditional  rectangular  grid.  The  range 
resolutions  (Ax,  Ay )  can  be  increased  in  regions  of  interest  without  increasing  the  expense 
of  computing  resources  for  the  case  of  curvilinear  nearly-orthogonal  grid2.  For  the 
traditional  rectangular  gridl,  the  horizontal  range  resolution  has  to  be  increased  from  121 
to  at  least  350  grid  points  in  order  to  obtain  similar  numerical  solutions  of  the  curvilinear 
grid2.  This  results  in  a  matrix  size  of  350  x  191  =  66,877  grid  points  vice  121  x  191  = 
23,1 1 1  grid  points  which  brings  the  amount  of  computing  expense  to  almost  three  times 
higher  than  the  computing  expense  of  the  curvilinear  grid2  for  a  similar  solution.  The  grid 
generation  technique  also  has  the  capability  of  producing  multi-block  grids  which  is  not  the 
scope  of  this  study.  However,  the  applications  and  relative  advantages  of  this  capability 
can  greatly  enhance  the  methods  of  solving  systems  of  fluid  flow  equations  in  coastal  ocean 
modeling  and  are  a  subject  for  further  exploration. 
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APPENDIX  A:  SALINITY  FIELD  SOLUTIONS 


Appendix  A  contains  the  simulations  of  salinity  field  numerical  solutions.  The 
simulations  are  numbered  from  Figure  A1  to  A21  and  are  being  referred  to  as  such  in  the 
text. 
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Figure  Al:  Gridl  -  Initial  Surface  Salinity 
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Figure  A2:  Grid2  -  Initial  Surface  Salinity 
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Figure  A4:  Grid!  -  200  days  Surface  Salinity 
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Figure  A5:  Grid2  -  200  days  Surface  Salinity 
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Figure  A13:  Grid!  -  200  days  Surface  Salinity 


days  Level  10  Salinity 


Figure  A 18:  Grid3  -  200  days  Level  10  Salinity 


Figure  A19:  Gridl  -  200  days  Level  15  Salinity 
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APPENDIX  B:  TEMPERATURE  FIELD  SOLUTIONS 


Appendix  B  contains  the  simulations  of  temperature  field  numerical  solutions.  The 
simulations  are  numbered  from  Figure  B1  to  B21  and  are  being  referred  to  as  such  in  the 
text. 
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Figure  Bl:  Grid!  -  Initial  Surface  Temperature 
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Figure  B2:  Grid2  -  Initial  Surface  Temperature 
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Figure  B3:  Grid3  -  Initial  Surface  Temperature 
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Figure  B4:  Grid!  -  200  days  Surface  Temperature 
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Figure  B6:  Grid3  -  200  days  Surface  Temperature 
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Figure  B7:  Gridl  -  Initial  Surface  Temperature 
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Figure  B8:  Grid2  -  Initial  Surface  Temperature 


Figure  B9:  Grid3  -  Initial  Surface  Temperature 
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Figure  BIO:  Grid!  - 100  days  Surface  Temperature 
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Figure  Bll:  Grid2  -  100  days  Surface  Temperature 


Figure  B12:  Grid3  - 100  days  Surface  Temperature 
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Figure  B13:  Gridl  -  200  days  Surface  Temperature 
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Figure  B15:  Grid3  -  200  days  Surface  Temperature 
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Figure  B18:  Grid3  -  200  days  Level  10  Temperature 


71 


Figure  B19:  Gridl  -  200  days  Level  15  Temperature 
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Figure  B21:  Grid3  -  200  days  Level  15  Temperature 
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